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Abstract 

Detonation of a tliree-dimensional reactive non-isotropic molecular crystal is modeled using 
molecular dynamics simulations. The detonation process is initiated by an impulse, followed by 
the creation of a stable fast reactive shock wave. The terminal shock velocity is independent of the 
initiation conditions. Further analysis shows supersonic propagation decoupled from the dynamics 
of the decomposed material left behind the shock front. The dependence of the shock velocity 
on crystal nonlinear compressibility resembles solitary behavior. These properties categorize the 
phenomena as a weak detonation. The dependence of the detonation wave on microscopic potential 
parameters was investigated. An increase in detonation velocity with the reaction exothermicity 
reaching a saturation value is observed. In all other respects the model crystal exhibits typical 
properties of a molecular crystal. 

PACS numbers: 62.50.Ef, 47.40.Rs, 82.40.Fp 
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I. INTRODUCTION 



Explosives are characterized by a detonation wave which propagates through the materiaL 
After initiation, the velocity of the detonation front reaches a steady state that exceeds the 
speed of sound in the material. The present paper is devoted to the analysis of a model solid 
explosive with the purpose of correlating the microscopic structure and interatomic forces to 
the bulk detonation properties. The investigation is based on a classical molecular dynamics 
simulation with a simple force field. The initial goal was to construct a first principle 
model which is able to qualitatively reproduce a stable detonation wave. The microscopic 
parameters considered are the crystal structure, the intermolecular forces that stabilize 
this structure, and the intramolecular potentials which yield the driving chemical reaction. 
The present investigation unravelled a new type of a solitary-like detonation wave which is 
directly driven by the one-step exothermic chemical decomposition. From a hydrodynamical 
perspective it can be classified as a weak detonation. 

A common theoretical framework for simulating detonations relies on a continuum picture 



of the material properties using a 
of mass, momentum and energy 



lydrodynamic description, and is based on the conservation 
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In the Chapman- Jouguet (CJ) model, the detonation wave is considered as a discontinuity 
between "unburnt" and "burnt" material, assuming that the transition from one state to 
another does occur instantaneously. It assumes an infinite reaction rate and requires only 
the knowledge of the thermodynamic equilibrium state of the burnt material to satisfy 
the conservation relations on the Hugoniot curve of compressed detonation products. The 
resulting detonation wave has zero-width reaction zone and should satisfy the minimum 
entropy solution of the Rankine-Hugoniot equations which corresponds to sonic condition for 
downstream flow (downstream velocity with respect to the front, which equals the isentropic 
sound speed in reacted material). Freely propagating one-dimensional detonation wave 
without the piston [unsupported detonation) attains the steady-state velocity satisfying the 
sonic condition. The CJ model is used extensively as an engineering tool enabling quick 
calculation of the detonation velocity and pressure in the design of explosive devices. 

An extension of the CJ model was proposed by Zeldovich, von Neumann, and Doring 
js-?] (ZND) incorporates finite-rate reaction kinetics resulting in a chemical reaction zone 
of finite length. It assumes that the unreacted explosive material behind the front is initially 
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compressed to a high pressure resulting in so-called von Neumann pressure spike. Then the 
chemical reactions are initiated behind the spike, leading to the exothermic energy release, 
temperature increase and expansion of the reaction products to the lower pressure. The 
structure and length of the reaction zone is determined by the equation of state of interme- 
diate products. After completion of the reactions, the products reach the final equilibrium 
state, at the end of the reaction zone. The final state determines the detonation velocity 
and depends only on the equation of state of the reaction products. For an unsupported 
detonation wave, this state satisfies the sonic condition as in the CJ model, and is called the 
CJ state. 

In the ZND model, besides the minimum-entropy (or, correspondingly, minimum-velocity) 
solution to the conservation equations, there are two other possible detonation propagation 
regimes. In a case of piston-supported detonation, the downstream flow can be subsonic in 
respect to the front and is called strong or overdriven detonation. Another type of possible 
hydrodynamical solutions is termed weak detonation. It is characterized by a downstream 
flow behind the detonation front which is supersonic in respect to the burnt material and 
effectively decouples from it (i.e. no disturbances from the burnt material can overtake the 
front). In addition, the weak detonation wave is characterized by only a moderate increase 
in density and pressure which is smaller than the CJ pressure on the detonation products 
Hugoniot. 

Zeldovich js] stresses that the weak detonation solution satisfies the boundary conditions 
of shock propagation and therefore is an admissible hydrodynamical solution. He argues 
that the reason that the weak detonation phenomena has not been observed experimentally 
is the unattainability of the steady state conditions of this solution. In particular, a typical 
ZND route from the unreacted shock-compressed state to the burnt material gets trapped 
in the stable CJ point of minimum entropy production. Von Neumann, however, has shown 
that if the Hugoniot curves of partially reacted products intersect one another then the 
weak detonation solution is possible {gJ. For example, explosives that have very rapid ini- 
tial exothermic decomposition followed by a slower endothermic reaction may exhibit such 
behavior sometimes called pathological detonation. It has been also speculated by Tarver, 
based on a hydrodynamical equations, that in porous materials weak detonation could be 
stable {9]. Though standard ZND model assumes that chemical reactions are triggered by 
the high temperature due to the initial shock compression, it is also possible to consider 
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an alternate ignition mechanism without prehminary shock heating at the von Neumann 



spike. Zeldovich 



10| discussed possible initiation of weak detonation by triggering chemical 



reactions with a sequence of sparks artificially fired along a tube. 

Can a first-principle molecular dynamics (MD) simulation converge to the hydrodynam- 
ical detonation model? 

There has been a continuous effort to develop MD simulation methods of explosives 
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12| . An interaction potential, known as the Reactive Empirical Bond Order potential 



(REBO), was introduced in Ref. [13[. This potential has a complex functional form and 
several parameters. Molecular dynamics simulations of detonation through a model 2D 
crystal using this potential, showed agreement with the predictions of the CJ model 14 1. 



A significant step was the development of a reactive force field which accounts for break- 



ing and forming chemical bonds during the passage of the shock wave through the solid 
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Simulations of actual explosives, such as RDX, PETN and TATP, have been attempted 
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23| . The main goal is to compare the simulation output to known experimental charac- 



teristics of the system, such as detonation velocity and final reaction products distributions. 
These studies were devoted to establish a realistic simulation scheme which converges to 
the framework of the hydrodynamic models. It was found that the goal of reaching steady 
state detonation conditions requires a major computational effort. The phenomenon stabi- 
lizes only in a mezoscopic scale (micrometers) which requires very large scale calculations 
including many millions of atoms. To establish such a method as a predictive tool, many 
such simulations should be performed and compared with experiment. A second round of 
improvements should then be applied to the force field. At present, this task is still in its 
infancy. 

Schemes to bridge the gap between the hydrodynamical description of detonation and the 



MD approach have been explored. The idea is to replace a group of mo^ 



ecules by a single 



mesoparticle with an internal thermodynamic degree of freedom [2J, |25|] or to describe a 
hydrodynamical cell by fictitious particles {26J. In both these schemes individual molecular 
properties are overlooked. 

In the present study a different MD approach was utilized. We limited the objective to 
establish a relationship between the forces governing the dynamics in the microscopic system 
and the macroscopic phenomena. The simple molecular force field employed is constructed 
to have only a small number of adjustable parameters. In a three dimensional model, 
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we observed stable detonation waves, their characteristics being independent of the initial 
conditions. It required a more thorough investigation to identify that these detonation waves 
are of a different character from those in the standard ZND model. 

A first indication of a new phenomena can be inferred from the one- dimensional (1-D) 
model of Peyrard et al. T^, constructed from a chain of unstable diatomic molecules. 
The dissociation reaction generated an accelerating detonation wave. In order to obtain 
a stable detonation wave artificial dissipation channels were added. Our analysis of their 
model revealed that a detonation wave in the 1-D chain resembles a solitary wave. Further 
understanding requires a full 3-D crystal model which should supply a natural dissipative 
mechanism. 

The initiation of detonation wave in solids is closely related to the propagation of a shock 
wave through a crystal lattice. Both waves are quasi-one-dimensional. Early MD simulations 



of a planar shock wave in a perfect fee crystal 28|, |29(| revealed the formation of a solitary- like 



train at the shock front over a wide range of shock velocities. It has oscillatory structure and 
exhibits significant deviation from the thermal equilibrium inside the train. In addition, the 
MD simulations of shocks in a perfect bcc crystal have shown 3^ that an isolated solitary 
wave can be emitted from the shock front and run ahead it at significantly faster speed 
than the shock front velocity. This links the present study to the subject of solitary waves 
propagation in a discrete lattice. Solitons are characterized by a group velocity which is 
proportional to the amplitude of the wave. For a one-dimensional lattice we should refer to 
the work of Toda who established a relationship between the microscopic parameter of the 



repulsive part of the potential and the group velocity [31|, |32| . Holian has performed MD 
simulations of shock waves in the 1-D Toda lattice to study the dependence of solitary train 
structure at the shock front on the interaction potential parameters j33i] . 

Can the energy release reaction occur directly at the shock front as in the CJ model? The 
prerequisite is a metastable molecule whose one-step exothermic decomposition is triggered 
by the shock front and its energy is immediately fed back into the detonation wave. By 
definition, shock wave propagation is faster than any linear elastic wave in the same material. 
As a result it will always confront fresh unperturbed molecules. For these molecules to 
contribute to the detonation event, their initial decomposition timescale should be similar 
or faster than the timescale determined by the front propagation velocity. Slower processes 
can take place behind the shock front. These processes may include thermal equilibration 
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and heat dissipation as well as slower chemical reactions which lead to the final product 
distribution. If the slower reactions are endothermic, the weak detonation wave can be 
initiated. 

In the present study we identified a diatomic molecule where the exothermic decompo- 
sition is sufficiently fast and can proceed directly at the detonation front. We observed a 
propagation of fast supersonic reactive wave initiated by the impulse impact in a three- 
dimensional perfect molecular crystal. We analyse the dependence of the detonation wave 
parameters on the interaction potential and the molecule decomposition kinetics. We find 
that the downstream fiow is supersonic, and the density increase behind the front is very 
small, which is typical for a weak detonation. Wc also observe the decoupling of supersonic 
front propagation from the dissipative dynamics of the decomposed material left behind the 
front. Our analysis suggests that the weak detonation wave in the simulated crystal resem- 
bles behavior of a solitary wave propagating at the supersonic speed through the lattice. 
The molecule decomposition is triggered at the detonation front where a significant thermal 
non-equilibrium exists. The released energy is directly channelled into acceleration of the 
atoms from the dissociated molecules pushing solitary detonation front forward. Remark- 
ably, the front propagation dynamics significantly depends on the relative orientation of the 
light and heavy atoms in the molecule with respect to the shock direction. 

The simulation results show a possibility of fast initiation of the molecule decomposition 
reactions at the supersonic solitary wave front in contrast to the standard ZND model 
where the reactions take place behind the shock front, initiated by the shock compression. 
Summarizing, the solitary front ignition can provide a mechanism for a direct transition to 
the weak detonation regime from the initial uncompressed state without going first to the 
high pressure state by the shock compression. 

II. THE COMPUTATIONAL MODEL 

An explosive is defined as a molecular crystal that can decompose to smaller particles, 
generating a stable detonation wave. The front of the detonation wave moves faster than 
the speed of the linear compression wave in the molecular crystal. Molecular dynamics 
simulation of detonations requires a scheme of molecular forces, initial geometry, boundary 
conditions and atomic masses. The model should be able to sustain a stable crystal structure. 
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Simulation of the system is based on the solution of the classical equations of motion with 
sufficient accuracy, using a modified version of the MD code MOLDY, which is described in 
appendix |X1 



A. The Reactive Molecule Model 



The model crystal is represented by a slab of diatomic molecules arranged in a crystal 
structure. Each of the atoms can represent also a group of atoms. The two effective par- 
ticles will be marked as N and C, and the masses are 47 amu for N and 15 for C. These 
notations and masses originated from nitromethane: the N corresponds to the NO2 group, 
and the C to the CH3 group. The difference between the masses is essential in such reactive 
molecular model j27|. However, no other similarity to nitromethane exists in our model. 
The chemical bond between N and C is designed to be metastable. When energized it can 
dissociate exothermally. The shape of the potential energy curve as a function of the N-C 
distance, is shown in figure [1] When the distance N-C is smaller than the barrier position 
the molecule is bound. However, when this distance is increased beyond this position the 
molecule disintegrates. In an energetic perspective, the molecule absorbs sufficient energy 
such that its internal energy is higher than the energy barrier. Separation of the N-C atoms 
results in decomposition of the molecule to its constituents accompanied by energy release. 
The evaluation of the potential and force during the trajectory calculations was carried out 
by generating and storing a table of the potential and force values at different N-C dis- 
tances, according to the potential shown in figure [H or similar potentials. A cubic spline 
interpolation was used to extract the values. There are some functional forms that can be 
used to generate such an exothermic reactive potential, some of them are described in Ref. 



We used a piecewise defined function to generate the potential, its functional forms are 



described in appendix iBl 



B. The Molecular Crystal 

The molecular crystal was constructed by a lattice of N-C molecules with a bonding 
interaction between similar groups in neighboring molecules. The N entity interacts with 
other N and C interacts with other C entities. A Morse potential was chosen to describe 
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FIG. 1: The reactive internal potential as a function of the N-C distance, for two different ratios 
between the barrier height and the exothermicity. When the distance N-C is smaller than the 
barrier position (here at 1.9 A - the local maximum) the molecule is bound. When the distance 
is increased beyond this position the molecule dissociates. Initially, the N-C distance is at the 
equilibrium distance (here at 1.5 A - the local minimum). Q represents the exothermicity, and Di 
the barrier height for the N-C dissociation. 



these interactions: 

V{r) = D (e-2^(^-"«) - 2e-^("-''«)) (1) 

The parameters used are shown in table [B The initial distances between nearest neighbor 
molecules was chosen as = 7. 5 A. A long and narrow slab of this molecular crystal 
was assembled with FCC symmetry. Other structures studied, the BCC and simple cubic 
crystals, were found to be unstable with this type of pair potentials. This finding is consistent 
with analysis of stability of Lennard- Jones crystals 34 1. 



C. Boundary Conditions and Initial Conditions 

We expect the detonation wave to be quasi one-dimensional. Hence, the reactive crystal 
slab used is chosen to be long and narrow. The direction of the propagation axis was chosen 
as Z, and it coincides with the [111] crystallographic axis. The molecular axes were also 
oriented along the Z direction. Thus, the crystal is very non-isotropic. The length of the 
periodic unit in the [111] direction of an FCC crystal is \/6 times tq. In our case tq = 7.5A, 
so the length of the periodic unit is ~ 18.37A. This distance will be referred to in the 
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D [eV] 


/3 [A-^] 


ro [A] 


N-N 


0.08 


1.0 


7.5 


C-C 


0.016 


1.0 


7.5 


N-C (inner part): 








barrier of 0.25 eV 


0.25 


4.33 


1.5 


barrier of 2.0 eV 


2.0 


4.33 


1.5 



TABLE I: Potential parameters used for the inter-molecular interaction. The third line shows the 
parameters of the inner part of the intramolecular reactive bonding potential (see appendix IB]) . 
The same potentials were used for N-C atoms of different molecules. 

following as the unit cell length. The length of the crystal was 295 unit cells along the Z 
direction, and, in some cases, a crystal with a length of 587 unit cells was used. Periodic 
boundary conditions were imposed along the perpendicular directions X and Y. To check 
that the size of the cross section is sufficient we ran simulations of crystals with various cross 
sections. It was found that detonation properties are converged in the larger cross sections. 
The converged cross section which was eventually used consists of 48 molecules in the XY 
plane. 

The shock wave was initiated by a small crystal pellet, assigned with an initial high 
velocity in the Z direction. The pellet collided with one of the small faces {XY plane) of 
the slab. The velocity of the pellet was not kept constant: the MD simulations were of 
NVE type, and the positions and velocities of the pellet's particles were calculated with no 
distinction. The pellet was composed of 3 crystal layers along the Z direction and had 48 N- 
C molecules in each layer. The collision of this pellet was sufficient to initiate a sustainable 
shock in the primary crystal. 

Most of the simulations were carried with initial temperature of 0°K. A few simulations 
were carried with higher initial temperatures. 

III. WEAK DETONATION WAVES IN THE CRYSTAL 

The analysis of weak detonation wave starts from the static properties of the reactive crys- 
tal. The second step describes the initiation process, showing that a stationary detonation 
wave is formed, independent of initiation process. The classification as a weak detonation 
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FIG. 2: A snapshot of the detonation wave. The propagation axis is the Z direction. Periodic 
boundary condition are used in the X, Y directions. The right hand side shows the unperturbed 
crystal structure. The red objects represent the heavy particles (N). The black objects represent 
the light particles (C). On the left side of the image one can observe the burnt material after the 
passage of the shock wave. The shock front is characterized by pilot cascades of light particles 
which are emitted from the decomposed molecules, and initiate the next layer in a domino-like 
effect. The lower panel is an enlarged viewpoint of the detonation front. The arrows indicate 
decomposed N-C pairs, corresponding to the pilot cascade. 

This simulation was carried out with exothermicity of 6 eV and a barrier of 0.25 eV. Eventually, 
in this simulation, in approximately 20 layers, most of the material is decomposed. 

relies on thermodynamics analysis of the phenomena. 

The equilibrium properties of the crystal model were characterized. The details are 
summarized in appendix O A linear relation was found between the velocity of the elastic 
p-wave and the stiffness of the potential (the /3 parameter in Eq. ([1]) ). 

A. Initiation of stable detonation vi^aves 

We ran NVE simulations of shock waves in the model reactive crystal, initiated by the 
small pellet described above. After initiating a shock wave, a transient emerges, eventually 
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stabilizing into a stable shock front. Such simulations were run with different crystal pa- 
rameters and different initial velocities of the pellet. The propagating shock wave initiated a 
decomposition reaction in crystals composed of reactive molecules. When the decomposition 
reaction kept pace with the shock front, the shock wave transformed into a detonation wave. 
The simulations were carried with initial temperature of 0°K. A few simulations were carried 
with higher initial temperatures, in order to verify that the phenomena is not restricted to 
0°K. 

A stable detonation wave is independent of initial conditions. This was verified by using 
different initial pellet velocities leading to shock waves that are independent of these initial 
velocities. The details can be found in appendix |Dl and in Fig. [T71 

A snapshot of the detonation process is shown in Fig. [2j A distinction between the 
unperturbed crystal and the burnt material is clearly seen. An enlarged section of the 
reaction zone is shown on the bottom. 

Similar simulations, with an endothermic intramolecular Morse potential, yielded decay- 
ing shock waves, with propagation velocity that depends on the initial pellet velocity. 



B. Classification of the Detonation Waves as Weak Detonations 

After the short transient the shock wave stabilized into a stable weak detonation wave. 
The shock front propagated in a very high Mach number, practically decoupled from other 
processes in the crystal. The number of molecules that were decomposed by this shock front 
changed from one simulation to another, depending on internal potential parameters such 
as the exothermicity and the barrier height for C-N dissociation. Only a small compression 
was associated with this front. Other properties, such as the mass velocity (Cf. IIIIC2l and 
the kinetic energy (Cf. IIII C 3|) . also exhibit an increase in magnitude. This shock wave will 



be referred to as the reaction wave in the following. 

The reaction wave travels through an unperturbed crystal, and therefore its velocity is 
determined by the crystal properties (cf. Sec. [IV|) . The reaction wave propagation velocity, 
as well as its amplitude, are independent of the initial impact (cf. appendix [D]). The 
velocity does depend on the exothermicity of the reaction, which is an intrinsic property of 
the material. The dependence of the detonation velocity on the reaction exothermicity will 
be discussed later (cf. IIV Al) . 
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After the reaction wave, another shock front appeared, characterized by a large increase 
of both the mass velocity and the kinetic energy. In the second shock front most of the 
remaining molecules dissociated, leading to a decomposition yield of above 85%. This second 
shock wave will be referred to as the compression wave in the following. In contrast to the 
reaction wave, the amplitude of the compression wave decays during its propagation. The 
compression wave velocity depends on the impact strength but not on the exothermicity. 

The reaction wave is faster than the compression wave: UreacUon > Ucompression (cf- Figs. 
OandH]). This fact categorizes the reaction wave as a weak detonation: The compression 
wave is a typical shock wave. It moves faster than sound waves in the reacted material, 
therefore it sets an upper limit for the sound wave velocity: Ucompression > Usound- Thus: 
Ureaction > Usound- Wc scc that the rcactiou wavc is supcrsouic with respect to the matter 
behind it. This is the hydrodynamical definition of weak detonation [4, p. 280]. 

C. Thermodynamics Profiles of the Weak Detonation Waves 

The profiles of weak detonations differ from those of normal detonations: In weak deto- 
nations there is no compression of the material before the reaction. After the reaction front, 
there is only a small change in the density, pressure, and the velocity of the particles. 

1. Different dissociation barrier height 

Fig. [3] compares the decomposition fraction of two simulations as a function of time and 
position. The first corresponds to molecules that have an N-C dissociation barrier of 2 eV 
(left plot). The second corresponds to molecules with a dissociation barrier height of 0.25 eV. 
The reaction exothermicity in both simulations is 6 eV. For a dissociation barrier of 2 eV, 
~ 35% of the molecules are decomposed by the reaction wave. When the compression wave 
passes through the partially "burnt" material, most of the remaining molecules decompose 
(small fluctuations around 90%). For a smaller dissociation barrier, 0.25 eV, more than 85% 
of the molecules are decomposed simultaneously with the reaction shock front. 

The reaction wave velocities are almost independent of the dissociation barrier height, as 
can be seen in Fig. [31 The barrier does influence the mass velocity and the decomposition 
fraction. The decoupling of the reaction wave from the compression wave and from the 
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FIG. 3: Contour plots of decomposition percentage during two different simulations. The left 
plot is obtained in a simulation with a dissociation barrier of 2 eV . Some of the molecules are 
decomposed by the reaction wave. Most of the other molecules (around 90%) are decomposed by 
the compression wave. The right plot is obtained in a simulation with a decreased dissociation 
barrier of 0.25 eV . More than 85% of the molecules are decomposed by the reaction wave. In both 
cases the exothermicity is 6 eV ^ and the initiating pellet velocity is ~ 90 km/ s. 
The reaction front velocity in the case of the higher barrier, ~ 60 km/ s, is very similar to the case 
of the smaller one, , ~ 62.5 km/ s. 



dissociation process results from the fact that this is the fastest wave in the material, leaving 
behind the slower processes. This phenomena is another signature of a weak detonation 
wave. 



2. Mass velocity profile 

The mass velocity characterizes the mass current, is defined as (v) = . Figure H] 

shows a contour plot of the Z component of the mass velocity during the simulation. The 
mass velocity plot shown can be compared to the decomposition percentage during the same 
simulation (right hand side of Fig. [3]): The decomposition front and the discontinuity of 
the mass velocity coincide, and define the reaction wave. The compression wave is observed 
on the mass velocity plot, but is absent from decomposition plot since the majority of the 
molecules have already decomposed. 

Figure [5] shows profiles of mass velocity at different snapshots during the simulation. The 
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FIG. 4: Contour plot of mass velocity along the Z direction. The dissociation barrier height 
in this simulation was 0.25 eV ^ and the reaction exothermicity was 6 eV . The initiating pellet 
velocity was ~ 90 km/s. Three regions are marked in the plot: Region A is the pre-shocked 
material. The mass velocity here is 0. Region B is the material after the reaction wave. The 
mass velocity here fluctuates around an average value of 0.3 km/s (see Fig [5] for more detailed 
profiles) . The mass velocity plot shown here can be compared to the decomposition percentage 
during the same simulation, which is shown on the right plot of figure [S) The decomposition front 
and the discontinuity of the mass velocity coincide, characterizing the reaction wave. Region C is 
the material after the compression wave. The mass velocity on the peak of the compression wave 
decays from around 20 km/s at 1 ps to 8.5 km/s at 10 ps. The compression wave is visible on the 
mass velocity plot, but is almost absent on the decomposition plot: For this barrier height, most 
of the molecules have already decomposed before the compression wave reached them. 

values were averaged over bins with a constant number of particles, 96, in each bin. The top 
plot shows the Z component of the mass velocity along the solid for three time values. The 
largest amplitude (the left peak for each time) is associated with the compression wave. At 
later times the peak mass velocity decays due to dissipation. The first shock front (the right 
sudden change in velocity for each time) is associated with the reaction wave. The insert 
shows a zoom- in on the reaction fronts. The average mass velocity after the reaction front, 
~ 0.3 km/s, is much smaller than the reaction front velocity, ~ 65 km/s. This is consistent 
with a small increase of density after the reaction front (Cf. Fig. [8]). It is clear that 
the reaction wave amplitude does not decay. When comparing simulations with different 
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FIG. 5: Mass velocity profiles along Z and X directions. The exothermicity of the reaction in 
this simulation was 6 eV ^ and the barrier was 0.25 eV . The values were determined by averaging 
over bins with constant 96 particles in each bin. The left edge of the crystal is placed at -500 A. 
The top plot shows the profile of the mass velocity in the Z direction, at three time values. The 
bottom plot presents the mass velocity along the X direction at a single time value. Along the Z 
direction (top figure), the largest amplitude (the left peak on each time step) is associated with 
the secondary compression wave. At later times the peak mass velocity decays due to dissipation. 
The first shock front (the right peak on each time step) is associated with the reaction wave. The 
reaction wave amplitude does not decay. The insert shows a zoom of the first shock fronts for the 
two shorter times. The average mass velocity after the reaction front, ~ 0.3 km/ is much smaller 
than the reaction front velocity, ~ 65 km/s. 
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parameters, we find that the reaction wave amplitude depends on the reaction exothermicity 
and is independent of the impact strength. As discussed above, in the introduction and in 
section [Dl this is a characteristic of a detonation wave. The compression wave behaves 
differently: it has a decaying amplitude that does not depend on the exothermicity but does 
depend on the impact strength. The compression wave is a regular compression shock wave. 
It travels through an unstructured matter composed of reaction decomposition products. 

The X component of the mass velocity is shown in the bottom plot for a single snapshot. 
The shock front is characterized with a abrupt onset of fluctuations. There is a minor change 
after the first shock front. 



\ 
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FIG. 6: Normalized kinetic energy profile at time = 10 ps. The exothermicity of the reaction in 
this simulation was 6 eV, and the dissociation barrier was 0.25 eV. The values were averaged over 
bins with constant number of particles, 96, in each bin. The blue line corresponds to kinetic energy 
along the X direction, and the red line to kinetic energy along the Z direction. The large peak on 
the right corresponds to the reaction wave front. 
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3. Energies and temperature 

Further insight on the detonation process can be obtained by inspection of the kinetic 
energy and the local temperature profiles, shown in figures |6] and [TJ respectively. The 
profiles correspond to an intermediate time of 10 ps. The local temperature is defined as 
Ti = ("^i ~ where {v) is the mass velocity defined above. The values were 

determined by averaging over bins with a constant (96) number of particles. The reaction 
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FIG. 7: Local temperature profile at time = 10 ps. The parameters are the same as in Fig[6l The 
blue line represents temperature along the X direction, and the red line along the Z direction. The 
insert shows an unthermalized region: The large peak on the right is the reaction wave front. The 
increase of the Z component of the local temperature precedes the increase of the X direction. 

wave front can clearly be seen as a sharp change in kinetic energy and temperature. The 
insert in the local temperature figure shows a small region after the first shock front which 
is not thermalized. The increase of the Z component of the local temperature precedes the 
corresponding increase of the X component. 

Figure |8] shows the potential energy and the density profiles at the same time step. On 
the right, a small increase of the potential energy is observed, followed by a large decrease. 
The small increase is caused by a minor compression before the dissociation starts, as seen 
in the density profile. The large decrease that follows is caused by the energy release of 
the exothermic decomposition process. The average density after the reaction front is only 
slightly larger than the density before this front. This is consistent with the small increase 
in the mass velocity (Cf. Fig. |5]). This phenomena characterizes weak detonations. 

IV. SATURATION OF THE WEAK DETONATION VELOCITY 

The detonation wave can be characterized by a stable detonation velocity independent 
of initiation parameters, as discussed in appendix |Dl The terminal velocity of the wave is a 
function of the microscopic parameters used in the model. Insight aimed at deciphering the 
phenomena is obtained by a systematic study of the variation of the detonation velocity as 
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FIG. 8: Potential energy (Ep) and density profiles at t = 10 ps. The parameters are the same as 
in Fig [6l The first shock front displays a small increase of the potential energy, due to a small 
compression before the dissociation. This small increase is followed by a large decrease of the 
potential energy, due to the energy release of the exothermic decomposition. The average density 
after the reaction front is not significant higher than the density before this front. The peak on 
the left is associated with the compression wave. 

a function of the parameters of interest. 



A. Intramolecular potential parameters: exothermicity of the reaction 

Simulations of detonation using reactive slabs were carried out with different values of 
the exothermicity. All other parameters, such as crystal structure, impact magnitude, in- 
termolecular potential, particles masses, initial conditions, were kept constant. Figure M 
displays the reaction propagation, which is determined by the dissociation front, in crys- 
tals with different exothermicity values. In the exothermicity range displayed in the figure 
(1.0 — 3.5 eV), the detonation velocity increases as a function of exothermicity. At larger 
exothermicity values, the detonation velocity reaches saturation: simulations with exother- 
micity values above 3.0 eV, show only a minor increase in the detonation velocity. Figure 
[TU] displays the variation in detonation velocity as a function of reaction exothermicity. The 
kinetic energy of the particles behind the reaction wave was found to depend linearly on the 
exothermicity of the reaction. Thus it is responsible for the energy balance. Nevertheless, in 
weak detonations this additional kinetic energy is decoupled from the reaction shock front. 
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FIG. 9: Reaction front propagation in crystals with different exothermicity values. All other 
parameters are identical. 
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FIG. 10: Detonation velocity in the crystal as a function of reaction exothermicity. The detonation 
velocity is measured in reduced units, normalized to the p-wave velocity, which is 4.6 km/ s. We can 
see the increase of detonation velocity as a function of exothermicity at low exothermicity values, 
until approximately 3 eV. The saturation of the detonation velocity at high exothermicity values 
can also be seen. The saturation was verified by simulations with very high exothermicity values: 
there is a very small increase of the detonation velocity in crystals with higher exothermicity values 
(12 eV, 15 eV, 18 eV, and even 27 eV. These results are not shown here). 



The saturation effect is an interplay between three factors: the crystal rigidity, responsible 
for the non-linear wave propagation, shear effects, responsible for energy dissipation at the 
reactive shock front, and kinematic effects, which determine the partitioning of the energy 
release during the decomposition of a single molecule. 
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Some insight regarding kinematic effects and the crystal rigidity is discussed in the fol- 
lowing, cf. sections IIVBI and IIV CI The non-linear effect of shear was not explored system- 
atically in this study. The convergence of the calculations was checked to an increase in the 
XY cross section (cf. Ill Cp . While checking convergence, we found that increasing the cross 
section decreases the detonation velocity. This is an indirect indication that the detonation 
velocity is sensitive to shear. 



B. Kinematic effects 



In weak detonations the decomposition mechanism of the reactive molecules has to be 
linked to the characteristics of the detonation wave. A propagation by a domino-like effect, 
as seen in Fig. [21 requires the decomposition to be asymmetric with respect to the Z 
direction. Insight is obtained by examining a simplified 1-D model of decomposition of a 
single diatomic molecule. Initially the molecule is at rest. A light particle emerges from the 
decomposition of the neighboring molecule, and hits the molecule from the heavier side (Fig. 
[TTl top). The collision initiates the decomposition process. At the end of the decomposition, 
the light particle of this molecule is emitted toward the next molecule (Fig. [TTl bottom). 
The mass of the light particle is denoted by m, and the mass of the heavier particle is am 
(a > 1). After the initiation collision, the molecule decomposes and the heat of the reaction, 

releases as kinetic energy. The velocity of the colliding particle is denoted by v\. After 
the decomposition, the velocities are denoted by Ui,U2,U3, for the colliding particle, the 
heavy particle, and the light particle, respectively. 

Under these conditions the equations for momentum and energy conservation become: 

mvi = mui + amu2 + mu^ 

^rnvf + E = ^rnvf + + \'mu\ 

where we assume that the collision is complete, meaning that all potential energy has been 
consumed. In steady state detonation, the initial velocity of the colliding particle should be 
equal to the velocity of the emitted light particle: 

vi = M3, (3) 

and this should also become the group detonation velocity. 
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FIG. 11: A molecule is hit by a light particle, emitted from neighboring decomposition (top 
diagram). After the collision, the molecule decomposes (bottom diagram). 

Two mechanisms come to mind: a direct and a delayed decomposition. In the direct 
decomposition mechanism, Eq. and give the complete kinematic picture. Substituting 
Eq. (E]) into Eq. ([2]) cancels out the terms Vi and U3. Therefore, the detonation velocity 
is undetermined without further assumptions. In order to determine the wave propagation 
velocity, more details on the dissipation mechanism, the crystal structure and the non-linear 
properties are needed. 

Next, we analyse the delayed decomposition mechanism. This mechanism has two steps: 
first, the colliding particle has an elastic encounter with the molecule. The molecule after the 
collision will acquire a velocity of U23 and mass of {a + l)m. On the second step, the molecule 
will dissociate into its components. The equations for momentum and energy conservation 
of the first step are: 

mvi = mui + (a + l)mu23 

(4) 

\mv\ = \mu\ + \{a + l)mu\.^ 
And the equations for momentum and energy conservation of the second step are: 

(a + l)mu23 = amu2 + mu?^ 

\{a + l)mu22, + E = ^amul + ^mul 

Again, we require the steady state condition, Eq. ([3]). Solving these equations shows that the 
velocities vi = M3, which are equal to the detonation velocity, are proportional to Ejm. 

A strong detonation process is adequately described by a delayed mechanism, since there 
is a delay between the compression (which leads to collisions) and the decomposition. The 
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yjE jm dependence of the strong detonation velocity is in accordance with the hydrodynam- 
ical theory. This a/E dependence was obtained in Ref. [ssl in MD simulations of detonations 



in the REBO model, and in Ref. 
tained in an earlier study 



in another crystal model. A different result was ob- 



371], in which a linear dependence of the detonation velocity on 



the energy released is reported. 

A weak detonation process is adequately described by the direct mechanism since the 
reaction coincides with the shock front. In weak detonations, we find the detonation velocity 
to saturate, becoming independent of the energy release. This means that in our 1-D model, 
only the results of a direct decomposition mechanism can fit this phenomena. The relation 
of direct reaction mechanism to weak detonations was shown also in hydrodynamic analysis 
of the ignition stage 



The mechanism discussed above suggests a difference between the two possible reaction 
wave propagation directions: Different behavior is expected if the reaction wave emits the 
heavier particle of the molecule, which in turn collides with the lighter particle of the neigh- 
boring molecule. This anisotropic behavior of the reaction wave was examined by initiating 
the shock wave in the opposite edge, impacting from the lighter side of the molecules. Under 
these conditions, no stable constant velocity reaction waves were formed. 

The dependence of the detonation wave velocity on the mass of the emitted particle 
was also examined. We performed simulations of detonations where the mass of the lighter 
particle was varied while maintaining the same structure and force field. In order to keep 
the total mass unchanged, the mass of the heavier particle was also changed. The masses 
used were 22, 10 and 5 amu for the C particle (and 40 ,52 and 57 amu for A^, respectively). 
Several simulations were performed, with different exothermicity values. We expect that 
decreasing the mass of the C lighter particle will result in a higher velocity, hence faster 
detonation. This has been already suggested by Peyrard et al. {2^ in a 1-D model. The 
magnitude of exothermicity at which the detonation velocity saturated was found in all cases 
to be close to 3 eV. However, the saturation velocity depended on the mass of the emitted C 
particle. Lower mass of this particle led to an increased detonation velocity. These results 
are presented in Fig. [T21 in log-log scale. A linear fit to the velocity logarithm as a function 
of the mass logarithm was calculated and is displayed in the graph. The resulting slope is 
close to 0.5 (0.49), which reveals a Xj ^fm dependence of the detonation velocity on the mass 
of the emitted group. This dependence is identical to the dependence of the shock velocity 
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of solitons in a Toda lattice. 
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FIG. 12: Saturation values of detonation velocities, marked with asterisks, as a function of the 
mass of the emitted C group. The graph is plotted in logarithmic scale on both axis. In all cases, 
this saturation value is obtained when the exothermicity is close to 3 eV . A linear fit of the data 
in the log-log scale is also plotted. The slope of the linear fit is close to —0.5 (—0.49), indicating 
that the detonation velocity is proportional to 1/y/m. 

C. Crystal rigidity: Inter molecular potential parameters 

Simulations of the propagation of small amplitude displacements in the crystal showed 
that the elastic P-wave velocity depends linearly on the /3 parameter of the intermolecular 
Morse potential (Eq. ([1])), i.e. the stiffness. This is in agreement with the propagation of 
elastic waves in a 1-D crystal model. See appendix IC 2bl and Fig. [16] for details. 

The shock wave velocity is crucially dependent on the rigidity of the crystal, which is 
determined by /3. To explore this dependence we carried out simulations of detonations in 
crystals with different /3 values. The first effect observed is that an increase in /3 leads to an 
increase in magnitude of the exothermicity for which the saturated reaction wave velocity 
is observed. The second effect found is that the saturation reaction wave velocity increased 
exponentially as a function of /3. Figure [13] shows the logarithm of the saturated reaction 
wave velocity as a function of /3. The exponential dependence can be clearly seen. The 
exponential dependence suggests that the saturation of the reaction wave velocity depends 
on the repulsive part of the intermolecular potential. This is different from the elastic p-wave 
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velocity which is governed by the shape of the potential near equilibrium and varies linearly 
function of /3. 
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FIG. 13: Comparison of reaction propagation velocities (scaled logarithmically) in 5 different 
crystals, /3 ranges from 1 A ^ to 1.4 A ^ (when /3 is the coefficient inside the Morse potential's 
exponent). The exponential dependence suggests that the saturation of the detonation velocity is 
caused by the repulsive part of the intermolecular potential. This is different from the elastic p- 
wave velocity dependence, which is governed by the behavior of the potential near the equilibrium 
point (see appendix IC 2 bl and Fig. [T6]) . 



To verify that the saturation of the reaction wave velocity is governed by the repulsive part 
of the intermolecular potential, we ran additional simulations. In these simulations we varied 
the stiffness of the inner part of the repulsive potential without altering other parameters 
of the crystal: The potential was combined from a repulsive part of Morse potential with 
13 = 1.2 A ^ and an attractive part of Morse potential with (3 = 1.0 A ^. The reaction shock 
wave velocities in these crystals were compared to the previous simulations. We found that 
the reaction waves velocities were almost identical to those with /3 = 1.2 A ^. To conclude 
the reactive shock- wave velocity is governed by the repulsive part, and independent of the 



soft long range part of the 



been suggested by Toda 



potential. This phenomena characterizes solitary waves, as has 
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3^ and by Rolfe et al. |39|. 
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V. DISCUSSION 



The purpose of this study was to bridge the gap between a first principle microscopic 
model of detonation and bulk hydrodynamical theories. The first task was to obtain a stable 
detonation wave independent of initial conditions. For this task we constructed a reactive 
molecular crystal model characterized by pair potentials. The equilibrium properties of 
the crystal are typical. It is stable at low temperatures and melts at temperatures which 
scale with the binding energy. The model crystal also possesses linear elastic waves. The 
detonation potential was added by making the molecule metastable to dissociation releasing 
a significant amount of energy. The model fulfilled the expectations and a stable detonation 
wave was identified. 

Further analysis, which compared the results obtained in the simulations to hydrodynam- 
ical theory, revealed a puzzling picture. The detonation wave did not have the characteristics 
of the common solutions of the ZND model. Some of the differences are: 

• The compression at the shock front was minor. 

• The chemical reaction coincided with the shock front. 

• The shock velocity was supersonic with respect to the burnt material left behind. 

Searching for a meeting point with hydrodynamical theory, we conclude that the phenom- 
ena we observed in the MD simulation is a weak detonation. Weak detonation is a possible 
hydrodynamical solution in which the shock wave is supersonic with respect to the material 
left behind. The characteristics of weak detonations are different from normal detonations: 
The weak detonation velocity is higher, and the pressure after completion of the reaction is 
smaller than in the normal detonation case. Also, in weak detonations there is no compres- 
sion of the material before the reaction. Zeldovich argues that such solutions are usually 
unattainable for substances that are initially inert Q]. However, he points out that weak 
detonations might occur if the chemical reaction would start in the initial state without pre- 

n 

liminary heating of the substance by the shock wave [40] . Dremin states that if self- ignition 



occurs at a pressure lower than the CJ pressure, a weak detonation wave is observed [41 1. 

Why are weak detonations attainable and stable in our case? The answer is in the 
kinematic behavior of the crystal. The shock propagates with characteristic of a nonlinear 
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solitary wave. The propagation velocity is determined by the repulsive part of the inter- 
atomic potential which is similar to solitons on a Toda lattice. The other kinematic property 
is the mass of the group that emitted from the dissociating molecule. When the shock front 
comes from the heavy side, the light particle is breaking out with the majority of the kinetic 
energy, initiating the reaction on the next crystal plane. This can be imagined as a shock 
propagation by the domino effect which never looks back. The ideal domino effect has no 
change in the density after the front. This is accompanied by a negligible increase in the 
mass velocity. The melting of the crystal and the equilibration of the burnt material lag be- 
hind the shock front and are decoupled from it due to the supersonic velocity. Remarkably, 
when the molecule orientation is the opposite (heavy particle is placed ahead of the light 
one in the shock direction) then the detonation does not reach a steady state. It explains 
the details of the solitary wave ignition mechanism. If the light particle is pushed forward 
at the high (supersonic) velocity after the decomposition, it quickly hits another molecule 
in the next crystal plane, allowing the reaction to supersonic speed due to the 

domino effect. On the contrary, if the heavy particle is placed ahead, it breaks out at much 
slower velocity, delaying the propagation of the decomposition reaction. 
These observations pose additional questions: 

• What are the conditions that are necessary to observe weak detonation in experiments? 

• What are the conditions that an MD simulation will reconstruct the standard solutions 
of the ZND model? 

It seems that a prerequisite for a stable weak detonation is a stiff molecular crystal which 
supports fast propagation of nonlinear shock waves. In addition, the crystal should be very 
non- isotropic. The molecules are oriented with the light particle pointing toward the prop- 
agation direction. We have carried out preliminary simulations with triatomic molecules 
with similar effects. Anisotropic shock initiation sensitivity has been observed in detona- 



tion of explosives composed o 
anisotropic sensitivity 
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single molecular crystals 



42|. Recent calculations reveal this 



44| . Weak detonations were shown to occur in the shock initi- 



ation process in the case of inhomogeneous explosives or inhomogeneous initial conditions. 
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47| . Weak detonations in mixtures which have nonmonotonic energy release were shown 



experimentally 



48|. It has been demonstrated that a quasi-steady form of weak detonation 



plays an integral role in describing shock-induced transition to detonation in an explosive 
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material 



49| . The transition to normal detonation was shown to occur effectively at the 



point where the weak detonation slows to the CJ velocity. In cases of very porous mate- 
rials it can happen that the decomposition wave will remain faster than the compression 
wave, stabilizing the weak detonation solution [q]. Finally, a weak detonation requires a 
decomposition reaction which can follow pace with the shock front even at relatively low 
temperatures. This dictates a time scale of a few tens of femtoseconds. 

The prerequisite for the MD simulation to generate the standard results of the Chapman- 
Jouguet or the more elaborate ZND model is a more complex and slower chemical reaction 
which can justify the quasi-equilibrium assumption. This influence of the speed and com- 



plexity of the reaction can be seen even in 1-D models: Elert et al. 



50| used three body 



interaction potentials in a 1-D model, anc 
cial frictional forces. The studies in Ref. 



K)t a stable detonation without introducing artifi- 



23| using realistic force fields aim at this direction. 



Nevertheless, the quasi-equilibrium assumption has been criticized [2l|. In the present study 
we find that most regions of the detonating crystal are well described by quasi-equilibrium 
assumptions, except the vicinity of the shock front. 
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Appendix A: MOLDY - The Molecular Dynamics Program 



The basic simulation program used to integrate the molecular dynamics equations was 



MOLDY 



5l| . MOLDY is suitable for the present purpose due to two primary reasons: 



The equations of motion are integrated in MOLDY using a modified version of the 



Beeman algorithm 
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52| . During detonation simulations there is rapid energy ex- 



change between potential and kinetic energy. Simulations carried out with the Verlet 
algorithm failed to conserve energy properly. The Beeman algorithm, with higher 
accuracy in velocities, was found to be adequate. 
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The neighbor hst in MOLDY is built using the hnked cell method (see, for example, 



53|). In shock wave simulations the system is not homogeneous: near the shock front 



there is a domain with high density, and the common neighbor list algorithm is not 



efficient for such a situation 



111- 



Two modifications of the MOLDY code were introduced to fit the requirements of deto- 
nation simulations: 

• MOLDY is constructed to calculate potential energy and forces from a wide set of 
common analytic potentials. We added an option to calculate the potentials and 
forces from a pre-prepared stored table, using cubic spline interpolation. 

• The initial velocities in MOLDY's simulations are sampled from the Maxwell Boltz- 
mann (MB) distribution. The pellet initial velocity were modified so that they could 
be preassigned, while the velocities of the slab atoms are sampled from the MB distri- 
bution. 



Appendix B: The functional form of the exothermic potential 

The interatomic potential described in figure [H can be constructed by a piecewise defined 
function, which is composed from three segments, each having a different functional form: 

Di [e-Mr-r^^n) -if + Q < r < r™„ 

V{r) = <j cgr^ + csr^ + cir + Cq r^i„ <r < r^ar (Bl) 

D2 (1 - (e-Z'^Cr-r.^.) _ i)2j ^^^^ ^ ^ ^ ^^^^ 

With rmin = 1-5 A, Tbar = 1.9 A and Vcut = 7 A refer to the positions of the local minimum, 
the barrier, and the cutoff, respectively. Di = 0.25 eV is the energy barrier height, Q is the 
energy release during decomposition of the molecule (the exothermicity of the reaction), and 
D2 should satisfy the requirement: D2 = Di + Q, so the function will be continuous. On the 
first segment there is a shifted Morse potential, so under small displacements, the molecule's 
behavior is similar to the behavior under standard Morse potential. The third segment is an 
inverted Morse potential, and it serves as the repulsive potential at the decomposed state. 
The role of the polynomial in the intermediate segment is to link between the two edge 
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segments. Therefore, the coefficients Cq, Ci, C2 and C3, are chosen to make the function and 
its first derivative continuous. 

We determined the polynomial coefficients by the requirements for continuity of the 
potential and its first derivative at the points r„iin and Thar (the derivative vanishes at these 
points). There are four requirements, so the coefficients are determined uniquely. Using Q = 
1.5 eV and Di — 0.25 eV (and, consequently, D2 — 1.75 eV), and using A as the length unit, 
we got C3 = -7.8125, C2 = 39.8437, Ci = -66.7969, and Cq = 38.4141, with the appropriate 
units. One can see that we can obtain continuity of the second derivative as well, by choosing 
^1 = y/3/{rbar - rmvn) = 4.3301 and ^2 = /SiV^i/£>2 = 1.6366 A"\ However, this 
is not essential, since the actual calculations during the simulations were done using cubic 
spline interpolation. This interpolation guarantees continuity of the second derivative. 

For other exothermicity values, (i.e. different Q), we used the same functional form on 
the first and second segments, with the same parameters' values, except Q and Cq, which 
determine the exothermicity. On the third segment, we determined the potential values in 
a way that maid only little gap of the function at Vcut- We also constructed potentials with 
different values for the barrier energy, using the same method. 

Appendix C: Crystal Characterization 

Before shock propagations are examined the equilibrium structure and characteristics has 
to be determined. This was carried out by evaluating the position correlation function at 
different temperatures. In addition, the acoustical sound velocity was calculated. 

1. Stability and Phase Transition Temperature 

Equilibrium NPT molecular dynamics simulations were performed under constant pres- 
sure of 1 bar, and at different temperatures. Under these conditions the stability of the N-C 
molecules as well as the stability of the crystal and its melting point were evaluated. 

We found that the molecules are stable in the temperature range between 0°K and 250°K. 
When the temperature was raised beyond 250°K, some of the N-C bonds started to dissoci- 
ate. This is consistent with the low dissociation barrier of the molecule under examination, 
which is 0.25 eV. 
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The radial distribution functions (RDF) of the N-N pairs were calculated. The results 
of the calculations for 50°K and 250°K are shown in figure dH The observed peaks are in 
positions reflecting the FCC lattice. The RDF verify that the crystal is solid up to 250°K. 
As the temperature is raised from 50°K to 250°K, The peaks of the RDF get broadened, 
and at 250°K the RDF does not vanish after the third peak, as a result of large fluctuations. 

On higher temperature, at 500°K, some of the molecules decompose and the slab melts. 
The RDF of the system at this temperature shows liquid behavior. In figure (TS] there is a 
plot of the RDF of N-N pairs on 500°K. There are no sharp peaks at distances greater than 
the distance of the first nearest neighbors peak. This character of the RDF is typical for the 
liquid phase. 
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FIG. 14: Radial distribution function of the crystal. Results are of a NPT simulations with 504 
molecules at pressure of 1 bar. The left plot corresponds to 50°K, and the right corresponds to 
250°K. 



2. Elastic Sound Velocity 

a. The elastic sound velocity of our crystal 

The MD simulation was used to determine velocity of the pressure P-wave in the crystal. 
A small longitudinal displacement of few crystal layers along the [111] direction was used 
to induce these waves. The perturbation propagation velocities was then determined. In 
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FIG. 15: Radial distribution function of a NPT simulations with 504 molecules at pressure of 1 
bar and temperature of 500°K. 

order to ensure that the displacements are small enough, the simulations were repeated with 
various displacement amplitudes. It was then verified that the propagation velocity is inde- 
pendent of the displacement amplitude. We found that the P-wave velocity is approximately 
2.6 unit cells per picosecond (^4.6x10'^ m/s). 



b. P-waves Velocity on Morse Crystals: A Discussion 

A simple 1-D model may give an insight regarding the P-wave velocity in crystals for 
which the pair interaction is described by the Morse function. In this model there is an 
infinitely long chain of masses. Let us consider first an infinite harmonic chain. In this 
case the interaction between neighboring masses is harmonic, with a spring constant a. We 
denote the distance between two neighboring masses /. The equation of motion for the n'th 
mass is: 

motn = a{Xn+l - 2Xn + Xn-l) (CI) 

If we assume that the time dependent position is given by: 

Xn{t) = (A+e+'*^"' + A.e-^'^"') cos(wt + 0) (C2) 
we get the dispersion relation: 

uj{k) = 2uJo sin C^) (C3) 



2 
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where we defined cuq = \foiJm. The p-wave velocity for this case is 



dk 



k=0 



1 l^^ 

WofCOS I — 



k=0 



Wq/ = \ —I 

m 



(C4) 



This model can be used to describe the harmonic approximation of small vibrations 
around the equilibrium distance in the Morse potential (and other potentials). If the vibra- 
tion amplitude is small enough for the harmonic approximation to hold, the propagation 
velocity of small perturbations should be independent of the perturbation amplitude. In such 
cases, we substitute a = 2/3^D, the second derivative near the minimum of the potential. 
Now the p-wave velocity is 

= .[^131 (C5) 
m \ m 

This result suggests that the p-wave velocity scales linearly with /3. 
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FIG. 16: P-wave velocity in the crystal as a function of the f3 parameter in the intermolecular 
Morse potential. Each velocity was calculated by simulating propagation of small displacements 
in a crystal, with the relevant f3 parameter, and marked in the figure by a cross. The dashed hne 
is an extrapolation to smaU /3 values assuming linear scaling. 

The P-wave velocity in this simple 1-D model is not directly comparable to the P-wave 
velocity in the 3-D fee slab, since the elastic waves on anisotropic crystal is determined 
by the elastic tensor a, which has several parameters 



54j |. However, some insights can be 



gained: The MD simulation was employed to evaluate the P-wave velocity along the [111] 
direction of our slab. The assumption of the velocity's linear scaling as a function of /3 is 
tested by simulating propagation of small displacements in different crystals, while changing 
(3. The results are shown in figure Uni marked by crosses. The dashed line in this figure is 
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FIG. 17: Reaction front location versus time in simulations with three different initial velocities 
of the pellet. The location was determined by a significant dissociation ratio. The results are 
grouped in three categories: The rightmost bundle represents five different simulations with the 
lowest pellet velocity used (wlO km/s). The central bundle corresponds to simulations (5 runs) 
with medium pellet velocity (w50 km/s), while the leftmost bundle corresponds to simulations 
with the largest pellet velocity (~120 km/s). (These high values of initial velocities are used for 
demonstration: pellet velocities of 10, 50, and 120 km/s give rise to initial transient shock waves 
that are, respectively, slower, similar, and higher than the developing stable detonation wave, which 
is approximately 65 km/s. Smaller values of pellet's velocities were used throughout most of this 
paper). The different simulations in each one of the three groups correspond to different initial 
conditions used (slightly different initial pellet-crystal distance). A short time after the pellet 
collision with the slab, different shock wave velocities develop in each simulation. The differences 
between the trajectories at short times are shown in the insert of the figure. A stable, steady state, 
velocity is reached after an additional propagation period. The shock velocity depends only on 
crystal parameters, which are identical in all simulations. 

an extrapolation of the P = 1 case, assuming linear scaling. We can see that our assumption 
is valid in a wide range of the f3 values. 
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Appendix D: The transition to stationary detonation waves 



A stable detonation wave is independent of initial conditions. To check this hypothesis 
different initial pellet velocities were tested. The location of the reaction front was defined 
as the first point where molecular decomposition is identified. Figure [17] shows the location 
of the reaction front as a function of time in simulations with three different initial pellet 
velocities. Inspection of these results shows that a short time after the pellet collision with 
the slab, different shock wave velocities develop in each simulation. This is a transient: a 
stable, steady state, velocity is reached after a short additional propagation period. The 
steady state shock velocity depends only on crystal parameters. Therefore, it is clear that the 
pressure wave induced by the pellet in this model crystal develops into a constant velocity 
detonation wave. A similar transition from initiation dependent velocity to a steady state 
velocity that depends only on crystal parameters was previously reported [55 1. 
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